library(INLA)
library(rgdal)
library(ggplot2)
library(viridis)
library(gridExtra)
library(raster)
library(ggspatial)
library(rgeos)

Data

# Xylella fastidiosa occurrence data in the demarcated area of Alicante (Spain)

data_Xf <- read.table("./Data/data_Xf.txt")

# Coordinate reference system

newproj <- "+proj=utm +zone=30 +ellps=GRS80 +units=m +no_defs"
## -- Boundaries -- ##

# Study area

bound <- readOGR('./Data/Boundaries', 'bound')

# Mountain barrier boundary

MB <- readOGR('./Data/Boundaries', 'bound_MB')

bound_MB <- gDifference(bound, MB, checkValidity = 2L)

# Continuous barrier boundary

CB <- readOGR('./Data/Boundaries', 'bound_CB')

bound_CB <- gDifference(bound, CB, checkValidity = 2L)

# Discontinuous barrier boundary

DB <- readOGR('./Data/Boundaries', 'bound_DB')

bound_DB <- gDifference(bound, DB, checkValidity = 2L)
## -- Remove data inside the barrier -- ##

# Mountain barrier

out.barrier = over(bound_MB,
                   SpatialPoints(cbind(data_Xf$X, data_Xf$Y), proj4string = CRS(newproj)),
                   returnList = T)[[1]]

data_M <- data_Xf[out.barrier,]

# write.table(data_M, './Data/data_Xf_M.txt')

# Perimeter barriers

out.barrier = over(bound_CB,
                   SpatialPoints(cbind(data_Xf$X, data_Xf$Y), proj4string = CRS(newproj)),
                   returnList = T)[[1]]

data_B <- data_Xf[out.barrier,]

# write.table(data_B, './Data/data_Xf_B.txt')
## -- Read data -- ##

# Mountain barrier data

data_M <- read.table("./Data/data_Xf_M.txt")

# Perimeter barriers data

data_B <- read.table("./Data/data_Xf_B.txt")

Models

## -- Mesh -- ##

# Stationary

mesh_ST <- inla.mesh.2d(
  boundary = bound,
  max.edge = c(750, 3500),
  offset = c(0.1,-0.1),
  cutoff   = 750
)

# Mountain barrier

mesh_MB <- inla.mesh.2d(
  boundary = bound_MB,
  max.edge = c(750, 3500),
  offset = c(0.1, -0.1),
  cutoff   = 750
)


# Continuous barrier

mesh_CB <- inla.mesh.2d(
  boundary = bound_CB,
  max.edge = c(750, 3500),
  offset = c(0.1, -0.1),
  cutoff   = 750
)

# Discontinuous barrier

mesh_DB <- inla.mesh.2d(
  boundary = bound_DB,
  max.edge = c(750, 3500),
  offset = c(0.1, -0.1),
  cutoff   = 750
)

1. Stationary model

data <- data_Xf
mesh <- mesh_ST
# SPDE model definition

size <- min(c(diff(range(data$X)), diff(range(data$Y))))
range0 <- size / 2

spde <- inla.spde2.pcmatern(
  # Mesh
  mesh = mesh,
  # P(practic.range < range0) = 0.5
  prior.range = c(range0, 0.5),
  # P(sigma > 10) = 0.01
  prior.sigma = c(10, 0.01)
)

# Projector matrix

A.est <- inla.spde.make.A(mesh, loc = cbind(data$X, data$Y))

# Data stack

stk.est <- inla.stack(
  data    = list(y = data$ResLab),
  A       = list(A.est, 1),
  effects = list(spatial = 1:mesh$n,
                 beta0 = rep(1, nrow(data))),
  tag     = 'est'
)
# Model fitting

formula.1 <- y ~ -1 + beta0 + f(spatial, model = spde)

model.est_ST <- inla(
  formula.1,
  data              = inla.stack.data(stk.est),
  family            = "binomial" ,
  control.compute   = list(
    dic              = TRUE,
    cpo              = TRUE,
    waic             = TRUE,
    return.marginals = TRUE
  ),
  control.predictor = list(A       = inla.stack.A(stk.est),
                           compute = TRUE),
  num.threads       = 2,
  verbose           = FALSE
)

saveRDS(model.est_ST, "./Models/ST_model_est.rds")
# Read saved model

model.est_ST <- readRDS("./Models/ST_model_est.rds")
# Projection on a grid

dxy <- apply(bbox(bound), 1, diff)
r <- dxy[1] / dxy[2]
m <- 120
proj.grid.mat <-
  inla.mesh.projector(
    mesh,
    xlim = bbox(bound)[1, ],
    ylim = bbox(bound)[2, ],
    dims = c(r, 1) * m
  )

# NA to the values outside boundary

ov <-
  over(
    SpatialPoints(proj.grid.mat$lattice$loc, bound@proj4string),
    SpatialPolygons(bound@polygons, proj4string = bound@proj4string)
  )

i.temp <- is.na(ov)

# Projector matrix to predict

A.pred <-
  inla.spde.make.A(mesh, loc = proj.grid.mat$lattice$loc[!i.temp, ])

# Stack to predict

stk.pred <- inla.stack(
  data      = list(y = NA),
  A       = list(A.pred, 1),
  effects = list(spatial = 1:mesh$n,
                 beta0 = rep(1, dim(A.pred)[1])),
  tag     = 'pred'
)

stk <- inla.stack(stk.est, stk.pred)
# Prediction

model.pred_ST <- inla(
  formula.1,
  data = inla.stack.data(stk),
  family = "binomial",
  control.predictor = list(
    A       = inla.stack.A(stk),
    compute = TRUE,
    link    = 1
  ),
  control.mode      = list(theta = model.est_ST$mode$theta,
                           restart = TRUE),
  num.threads       = 3,
  verbose  = TRUE
)

saveRDS(model.pred_ST, "./Models/ST_model_pred.rds")
# Read saved model

model.pred_ST <- readRDS("./Models/ST_model_pred.rds")
# Index of the predicted values

idx <- inla.stack.index(stk, 'pred')$data

# Matrix to visualize

prob.mean_ST <-
  prob.sd_ST <- matrix(NA, proj.grid.mat$lattice$dims[1],
                       proj.grid.mat$lattice$dims[2])
prob.mean_ST[!i.temp] <-
  c(model.pred_ST$summary.fitted.values$mean[idx])
prob.sd_ST[!i.temp] <-
  c(model.pred_ST$summary.fitted.values$sd[idx])

Results

  • Posterior marginal distribution of parameters
summary(model.est_ST)$fixed
##        mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## beta0 -1.68 0.245     -2.205   -1.666     -1.233 -1.641   0

 

  • Posterior marginal distribution of hyperparameters
summary(model.est_ST)$hyperpar

 

  • Mean and SD of the posterior predictive distribution

 

# Raster of the posterior predictive mean for the difference between the models

mean_pred_ST <- matrix_to_raster(prob.mean_ST, proj.grid.mat = proj.grid.mat)
mean_pred_ST<-mask(mean_pred_ST, bound)

2. Mountain barrier model

data <- data_M
bound_barriers <- bound_MB
mesh <- mesh_MB

# Polygon with the barrier

in.tri <- inla.over_sp_mesh(bound_barriers,
                            y = mesh,
                            type = "centroid",
                            ignore.CRS = TRUE)
num.tri <- length(mesh$graph$tv[, 1])
barrier.tri <- setdiff(1:num.tri, in.tri)
poly.barrier <- inla.barrier.polygon(mesh,
                                     barrier.triangles = barrier.tri)

# SPDE model definition

size <- min(c(diff(range(data$X)), diff(range(data$Y))))
range0 <- size / 2

barrier.model <- inla.barrier.pcmatern(
  mesh,
  barrier.triangles =
    barrier.tri,
  prior.range = c(range0, 0.5),
  prior.sigma = c(10, 0.01)
)


# Projector matrix

A.est <- inla.spde.make.A(mesh, loc = cbind(data$X, data$Y))

# Data stack

stk.est <- inla.stack(
  data    = list(y = data$ResLab),
  A       = list(A.est, 1),
  effects = list(spatial = 1:mesh$n,
                 beta0 = rep(1, nrow(data))),
  tag     = 'est'
)
# Model fitting

formula.1 <- y ~ -1 + beta0 + f(spatial, model = barrier.model)

model.est_MB <- inla(
  formula.1,
  data              = inla.stack.data(stk.est),
  family            = "binomial" ,
  control.compute   = list(
    dic              = TRUE,
    cpo              = TRUE,
    waic             = TRUE,
    return.marginals = TRUE
  ),
  control.predictor = list(A       = inla.stack.A(stk.est),
                           compute = TRUE),
  num.threads       = 2,
  verbose           = FALSE
)


saveRDS(model.est_MB, "./Models/MB_model_est.rds")
# Read saved model

model.est_MB <- readRDS("./Models/MB_model_est.rds")
# Projection on a grid

dxy <- apply(bbox(bound_barriers), 1, diff)
r <- dxy[1] / dxy[2]
m <- 120
proj.grid.mat <-
  inla.mesh.projector(
    mesh,
    xlim = bbox(bound_barriers)[1, ],
    ylim = bbox(bound_barriers)[2, ],
    dims = c(r, 1) * m
  )

# NA to the values outside boundary

ov <-
  over(
    SpatialPoints(proj.grid.mat$lattice$loc, bound@proj4string),
    SpatialPolygons(bound@polygons, proj4string = bound@proj4string)
  )

i.temp <- is.na(ov)

# Projector matrix to predict

A.pred <-
  inla.spde.make.A(mesh, loc = proj.grid.mat$lattice$loc[!i.temp, ])

# Stack to predict

stk.pred <- inla.stack(
  data      = list(y = NA),
  A       = list(A.pred, 1),
  effects = list(spatial = 1:mesh$n,
                 beta0 = rep(1, dim(A.pred)[1])),
  tag     = 'pred'
)

stk <- inla.stack(stk.est, stk.pred)
# Prediction

model.pred_MB <- inla(
  formula.1,
  data = inla.stack.data(stk),
  family = "binomial",
  control.predictor = list(
    A       = inla.stack.A(stk),
    compute = TRUE,
    link    = 1
  ),
  control.mode      = list(theta = model.est_MB$mode$theta,
                           restart = TRUE),
  num.threads       = 3,
  verbose  = TRUE
)

saveRDS(model.pred_MB, "./Models/MB_model_pred.rds")
# Read saved model

model.pred_MB <- readRDS("./Models/MB_model_pred.rds")
# Index of the predicted values

idx <- inla.stack.index(stk, 'pred')$data

# Matrix to visualize

prob.mean_MB <-
  prob.sd_MB <- matrix(NA, proj.grid.mat$lattice$dims[1],
                       proj.grid.mat$lattice$dims[2])
prob.mean_MB[!i.temp] <-
  c(model.pred_MB$summary.fitted.values$mean[idx])
prob.sd_MB[!i.temp] <-
  c(model.pred_MB$summary.fitted.values$sd[idx])

Results

  • Posterior marginal distribution of parameters
summary(model.est_MB)$fixed
##         mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## beta0 -1.612 0.227     -2.093   -1.601     -1.193 -1.581   0

 

  • Posterior marginal distribution of hyperparameters
summary(model.est_MB)$hyperpar

 

Standard deviation of the spatial effect \(= e^{\theta_1}\)

Range \(= e^{\theta_2}\)

pander::pander(exp(model.est_MB$summary.hyperpar[,c(1,3,4,5)]), row.names=c("SD", "Range"))
  mean 0.025quant 0.5quant 0.975quant
SD 1.433 1.205 1.432 1.71
Range 3861 2919 3844 5212

 

  • Mean and SD of the posterior predictive distribution

 

# Raster of the posterior predictive mean for the difference between the models.

mean_pred_MB <-
  matrix_to_raster(prob.mean_MB, proj.grid.mat = proj.grid.mat)
mean_pred_MB <- mask(mean_pred_MB, bound)

3. Continuous barrier model

data <- data_B
bound_barriers <- bound_CB
mesh <- mesh_CB

# Polygon with the barrier

in.tri <- inla.over_sp_mesh(bound_barriers,
                            y = mesh,
                            type = "centroid",
                            ignore.CRS = TRUE)
num.tri <- length(mesh$graph$tv[, 1])
barrier.tri <- setdiff(1:num.tri, in.tri)
poly.barrier <- inla.barrier.polygon(mesh,
                                     barrier.triangles = barrier.tri)

# SPDE model definition

size <- min(c(diff(range(data$X)), diff(range(data$Y))))
range0 <- size / 2

barrier.model <- inla.barrier.pcmatern(
  mesh,
  barrier.triangles =
    barrier.tri,
  prior.range = c(range0, 0.5),
  prior.sigma = c(10, 0.01)
)


# Projector matrix

A.est <- inla.spde.make.A(mesh, loc = cbind(data$X, data$Y))

# Data stack

stk.est <- inla.stack(
  data    = list(y = data$ResLab),
  A       = list(A.est, 1),
  effects = list(spatial = 1:mesh$n,
                 beta0 = rep(1, nrow(data))),
  tag     = 'est'
)
# Model fitting

formula.1 <- y ~ -1 + beta0 + f(spatial, model = barrier.model)

model.est_CB <- inla(
  formula.1,
  data              = inla.stack.data(stk.est),
  family            = "binomial" ,
  control.compute   = list(
    dic              = TRUE,
    cpo              = TRUE,
    waic             = TRUE,
    return.marginals = TRUE
  ),
  control.predictor = list(A       = inla.stack.A(stk.est),
                           compute = TRUE),
  num.threads       = 2,
  verbose           = FALSE
)


saveRDS(model.est_CB, "./Models/CB_model_est.rds")
# Read saved model

model.est_CB <- readRDS("./Models/CB_model_est.rds")
# Projection on a grid

dxy <- apply(bbox(bound_barriers), 1, diff)
r <- dxy[1] / dxy[2]
m <- 120
proj.grid.mat <-
  inla.mesh.projector(
    mesh,
    xlim = bbox(bound_barriers)[1, ],
    ylim = bbox(bound_barriers)[2, ],
    dims = c(r, 1) * m
  )

# NA to the values outside boundary

ov <-
  over(
    SpatialPoints(proj.grid.mat$lattice$loc, bound@proj4string),
    SpatialPolygons(bound@polygons, proj4string = bound@proj4string)
  )

i.temp <- is.na(ov)

# Projector matrix to predict

A.pred <-
  inla.spde.make.A(mesh, loc = proj.grid.mat$lattice$loc[!i.temp, ])

# Stack to predict

stk.pred <- inla.stack(
  data      = list(y = NA),
  A       = list(A.pred, 1),
  effects = list(spatial = 1:mesh$n,
                 beta0 = rep(1, dim(A.pred)[1])),
  tag     = 'pred'
)

stk <- inla.stack(stk.est, stk.pred)
# Prediction

model.pred_CB <- inla(
  formula.1,
  data = inla.stack.data(stk),
  family = "binomial",
  control.predictor = list(
    A       = inla.stack.A(stk),
    compute = TRUE,
    link    = 1
  ),
  control.mode      = list(theta = model.est_CB$mode$theta,
                           restart = TRUE),
  num.threads       = 3,
  verbose  = TRUE
)

saveRDS(model.pred_CB, "./Models/CB_model_pred.rds")
# Read saved model

model.pred_CB <- readRDS("./Models/CB_model_pred.rds")
# Index of the predicted values

idx <- inla.stack.index(stk, 'pred')$data

# Matrix to visualize

prob.mean_CB <-
  prob.sd_CB <- matrix(NA, proj.grid.mat$lattice$dims[1],
                       proj.grid.mat$lattice$dims[2])
prob.mean_CB[!i.temp] <-
  c(model.pred_CB$summary.fitted.values$mean[idx])
prob.sd_CB[!i.temp] <-
  c(model.pred_CB$summary.fitted.values$sd[idx])

Results

  • Posterior marginal distribution of parameters
summary(model.est_CB)$fixed
##        mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## beta0 -1.79 0.375     -2.633   -1.755     -1.149 -1.691   0

 

  • Posterior marginal distribution of hyperparameters
summary(model.est_CB)$hyperpar

 

Standard deviation of the spatial effect \(= e^{\theta_1}\)

Range \(= e^{\theta_2}\)

pander::pander(exp(model.est_CB$summary.hyperpar[,c(1,3,4,5)]), row.names=c("SD", "Range"))
  mean 0.025quant 0.5quant 0.975quant
SD 1.496 1.202 1.494 1.875
Range 6141 4296 6103 9043

 

  • Mean and SD of the posterior predictive distribution

# Raster of the posterior predictive mean for the difference between the models.

mean_pred_CB <-
  matrix_to_raster(prob.mean_CB, proj.grid.mat = proj.grid.mat)
mean_pred_CB <- mask(mean_pred_CB, bound)

4. Discontinuous barrier model

data <- data_B
bound_barriers <- bound_DB
mesh <- mesh_DB

# Polygon with the barrier

in.tri <- inla.over_sp_mesh(bound_barriers,
                            y = mesh,
                            type = "centroid",
                            ignore.CRS = TRUE)
num.tri <- length(mesh$graph$tv[, 1])
barrier.tri <- setdiff(1:num.tri, in.tri)
poly.barrier <- inla.barrier.polygon(mesh,
                                     barrier.triangles = barrier.tri)

# SPDE model definition

size <- min(c(diff(range(data$X)), diff(range(data$Y))))
range0 <- size / 2

barrier.model <- inla.barrier.pcmatern(
  mesh,
  barrier.triangles =
    barrier.tri,
  prior.range = c(range0, 0.5),
  prior.sigma = c(10, 0.01)
)


# Projector matrix

A.est <- inla.spde.make.A(mesh, loc = cbind(data$X, data$Y))

# Data stack

stk.est <- inla.stack(
  data    = list(y = data$ResLab),
  A       = list(A.est, 1),
  effects = list(spatial = 1:mesh$n,
                 beta0 = rep(1, nrow(data))),
  tag     = 'est'
)
# Model fitting

formula.1 <- y ~ -1 + beta0 + f(spatial, model = barrier.model)

model.est_DB <- inla(
  formula.1,
  data              = inla.stack.data(stk.est),
  family            = "binomial" ,
  control.compute   = list(
    dic              = TRUE,
    cpo              = TRUE,
    waic             = TRUE,
    return.marginals = TRUE
  ),
  control.predictor = list(A       = inla.stack.A(stk.est),
                           compute = TRUE),
  num.threads       = 2,
  verbose           = FALSE
)


saveRDS(model.est_DB, "./Models/DB_model_est.rds")
# Read saved model

model.est_DB <- readRDS("./Models/DB_model_est.rds")
# Projection on a grid

dxy <- apply(bbox(bound_barriers), 1, diff)
r <- dxy[1] / dxy[2]
m <- 120
proj.grid.mat <-
  inla.mesh.projector(
    mesh,
    xlim = bbox(bound_barriers)[1, ],
    ylim = bbox(bound_barriers)[2, ],
    dims = c(r, 1) * m
  )

# NA to the values outside boundary

ov <-
  over(
    SpatialPoints(proj.grid.mat$lattice$loc, bound@proj4string),
    SpatialPolygons(bound@polygons, proj4string = bound@proj4string)
  )

i.temp <- is.na(ov)

# Projector matrix to predict

A.pred <-
  inla.spde.make.A(mesh, loc = proj.grid.mat$lattice$loc[!i.temp, ])

# Stack to predict

stk.pred <- inla.stack(
  data      = list(y = NA),
  A       = list(A.pred, 1),
  effects = list(spatial = 1:mesh$n,
                 beta0 = rep(1, dim(A.pred)[1])),
  tag     = 'pred'
)

stk <- inla.stack(stk.est, stk.pred)
# Prediction

model.pred_DB <- inla(
  formula.1,
  data = inla.stack.data(stk),
  family = "binomial",
  control.predictor = list(
    A       = inla.stack.A(stk),
    compute = TRUE,
    link    = 1
  ),
  control.mode      = list(theta = model.est_DB$mode$theta,
                           restart = TRUE),
  num.threads       = 3,
  verbose  = TRUE
)

saveRDS(model.pred_DB, "./Models/DB_model_pred.rds")
# Read saved model

model.pred_DB <- readRDS("./Models/DB_model_pred.rds")
# Index of the predicted values

idx <- inla.stack.index(stk, 'pred')$data

# Matrix to visualize

prob.mean_DB <-
  prob.sd_DB <- matrix(NA, proj.grid.mat$lattice$dims[1],
                       proj.grid.mat$lattice$dims[2])
prob.mean_DB[!i.temp] <-
  c(model.pred_DB$summary.fitted.values$mean[idx])
prob.sd_DB[!i.temp] <-
  c(model.pred_DB$summary.fitted.values$sd[idx])

Results

  • Posterior marginal distribution of parameters
summary(model.est_CB)$fixed
##        mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## beta0 -1.79 0.375     -2.633   -1.755     -1.149 -1.691   0

 

  • Posterior marginal distribution of hyperparameters
summary(model.est_CB)$hyperpar

 

Standard deviation of the spatial effect \(= e^{\theta_1}\)

Range \(= e^{\theta_2}\)

pander::pander(exp(model.est_CB$summary.hyperpar[,c(1,3,4,5)]), row.names=c("SD", "Range"))
  mean 0.025quant 0.5quant 0.975quant
SD 1.496 1.202 1.494 1.875
Range 6141 4296 6103 9043

 

  • Mean and SD of the posterior predictive distribution

# Raster of the posterior predictive mean for the difference between the models.

mean_pred_DB <-
  matrix_to_raster(prob.mean_DB, proj.grid.mat = proj.grid.mat)
mean_pred_DB <- mask(mean_pred_DB, bound)

Differences of the posterior predictive mean